Transcriptomic Analysis of MSTN Knockout in the Early Differentiation of Chicken Fetal Myoblasts

In mammals, Myostatin (MSTN) is a known negative regulator of muscle growth and development, but its role in birds is poorly understood. To investigate the molecular mechanism of MSTN on muscle growth and development in chickens, we knocked out MSTN in chicken fetal myoblasts (CFMs) and sequenced the mRNA transcriptomes. The amplicon sequencing results show that the editing efficiency of the cells was 76%. The transcriptomic results showed that 296 differentially expressed genes were generated after down-regulation of MSTN, including angiotensin I converting enzyme (ACE), extracellular fatty acid-binding protein (EXFABP) and troponin T1, slow skeletal type (TNNT1). These genes are closely associated with myoblast differentiation, muscle growth and energy metabolism. Subsequent enrichment analysis showed that DEGs of CFMs were related to MAPK, Pl3K/Akt, and STAT3 signaling pathways. The MAPK and Pl3K/Akt signaling pathways are two of the three known signaling pathways involved in the biological effects of MSTN in mammals, and the STAT3 pathway is also significantly enriched in MSTN knock out chicken leg muscles. The results of this study will help to understand the possible molecular mechanism of MSTN regulating the early differentiation of CFMs and lay a foundation for further research on the molecular mechanism of MSTN involvement in muscle growth and development.


Introduction
The chicken has been used as a model research organism for many years, and it also provides large amounts of protein for the human diet [1]. In the domestic poultry industry, increasing meat production is an important goal; therefore, a full understanding of the regulatory mechanisms of muscle proliferation and differentiation is needed for improvement of growth and meat production traits [2]. However, knowledge is limited regarding the molecular mechanisms underlying muscle growth and development in birds in general, and in broiler chickens in particular.
Myostatin (MSTN), also known as a growth factor and differentiation factor-8 (GDF-8), is mainly expressed in skeletal muscle cells [3], and it acts as an essential negative regulator of skeletal muscle growth in livestock [4]. MSTN was first discovered and confirmed as a negative regulator of muscle growth in mice [3], and its inhibitory effect was subsequently further confirmed in cattle, pigs, chickens, goats, and humans [5][6][7][8][9][10]. Previous studies in mice and other mammals have shown that MSTN exerts its biological effects through multiple signaling pathways, including the mitogen-activated protein kinase (MAPK) signaling pathway [11], TGF-β [12], and the PI3K/Akt [13] signaling pathway. Unlike the case for mammals, however, the functional and molecular mechanisms of MSTN in chicken muscle have not been well studied.
The limited research on poultry has mainly focused on statistical genomic methods that have examined correlations between genomic SNPs or mutations in the MSTN gene and growth traits [14,15]. Extending these types of studies to include functional genomics has been difficult because gene knockout in chickens or other birds is difficult to perform due to differences in the avian reproductive system and embryonic development [16]. Therefore, few experimental results related to the developmental effects of MSTN are available for chickens.
The development of skeletal muscle is determined by myogenesis, and the number of muscle fibers formed mainly depends on the degree of myoblast proliferation and differentiation [17]. Researchers have used CRISPR/Cas9 technology to knock out the MSTN gene in mouse myoblasts (C2C12 cells) and performed RNA-Seq and miRNA-Seq to explore the molecular mechanism of MSTN regulation of muscle cell proliferation [18]. Other studies have used CRISPR/Cas9 technology to generate MSTN knockout (KO) avian (quail) myoblasts (QM7) to examine the MSTN regulatory networks through RNA-Seq [19].
In the current study, we used the CRISPR/Cas9 system to knock out the MSTN gene in chicken fetal myoblasts (CFMs), which-like C2C12 and QM7 cells-are an important model for studying the molecular mechanisms of muscle growth and development [20,21]. Then we performed RNA-Seq technology on knockout (KO) and wild-type (WT) CFMs and used PCA analysis (a mathematical algorithm to reduce data complexity) to evaluate the similarities and differences between samples and determine that the samples can be grouped [22]. Finally, we screened MSTN-related genes, pathways, and biological processes in knockout (KO) vs. wild-type (WT) CFMs. Our research provides valuable evidence for the possible molecular mechanism by which MSTN regulates the early differentiation of CFMs and lays a foundation for further research on the molecular mechanism of MSTN involvement in muscle growth and development.

Cell Culture
Specific pathogen-free (SPF) chicken eggs were obtained from Sais Poultry Co. Ltd., Jinan, China, and incubated at 37.8 • C under a relative humidity of 60%.
Chicken fetal myoblasts (CFMs) were isolated from chicken embryos (n = 20) at 13 days of age and minced in 0.2% collagenase type I (Solarbio, Beijing, China) for digestion. The resulting suspension was dispersed, filtered, centrifuged, and plated in cell culture plates. The cultures were enriched in myoblasts by placing the cells in an incubator twice for 40 minutes to remove fibroblasts. The enriched culture was grown in DMEM/F12 medium (Gibco, Shanghai, China) supplemented with 20% FBS (Gibco, Shanghai, China) at 37 • C 5% CO 2 . The CFMs were induced to differentiate by replacing the growth medium containing 20% fetal bovine serum with differentiation medium containing 2% horse serum.
As non-muscle control cells, primary chicken embryonic fibroblasts (CEFs) were isolated from chicken embryos (n = 8) at 9 days of age. The heads, wings, legs, and abdominal organs were removed from the chicken embryos, and the remaining tissues were minced and digested with 0.25% trypsin (Gibco, Shanghai, China). The resulting suspension was filtered, centrifuged, and plated in cell culture plates. CEFs were grown in DMEM/F12 medium (Gibco, Shanghai, China) supplemented with 5% FBS (Gibco, Shanghai, China) at 37 • C 5% CO 2 .

Immunofluorescence
Only primary myoblasts were used for immunofluorescence identification. Primary myoblasts were fixed with 4% formaldehyde (Solarbio, Beijing, China) for 30 min and then washed three times in phosphate-buffered saline (PBS) (Gibco, Shanghai, China). The fixed cells were then permeabilized with 0.2% Triton X-100 (Solarbio, Beijing, China) for 15 min and then washed three times in PBS. The cells were blocked with goat serum (Solarbio, Beijing, China) for 30 min and incubated with anti-desmin (Bioss, Beijing, China) overnight at 4 • C. The cells were then incubated with goat anti-rabbit IgG/FITC (Bioss, Beijing, China) for 1 h at 37 • C. The cell nuclei were visualized using DAPI staining solution.

Cell Transduction
Gently drop the AdV-EGFP system and the AdV-EGFP-CRISPR system into a six-well plate containing chicken CFMs, with a confluence of 50-70%. The cells were transduced at a multiplicity of infection (MOI) of 1000, with the addition of 5 µg/mL polybrene. After incubation for 12h, the cells were gently washed twice with PBS, and fresh growth medium was added. The same experiment was also performed on CEFs. One day after transduction, the growth medium of CFMs was replaced with differentiation medium so that the primary myoblasts began to differentiate. The CEFs were again cultured in growth medium. Two days after transduction, the cells were evaluated by fluorescence microscopy under appropriate excitation filters and cell samples were collected for subsequent analysis.

Detection of Editing Effect and Efficiency by Amplicon Sequencing
Two days after transduction, genomic DNA was isolated for PCR using the following primers: Primer F: CTGGATGGCAGTAGTCAGCC; and Primer R: CTGTTGGGAGAGC-CTGAGAA. PCR conditions were: 95 • C for 3 min; 35× (95 • C for 15 s, 56 • C for 15 s, 72 • C for 36 s); 72 • C for 5 min; followed by 4 • C. PCR products were shown on a 1% agarose gel. Experiment-specific barcodes were added to the 5' end of the primer sequence (Table S1) and performed PCR amplification. The resulting PCR products were pooled and sequenced with 250 bp paired end reads on a NovaSeq instrument. The samples were demultiplexed according to the assigned barcode sequences. Data were analyzed using the CRISPResso2 software. The alignment of reads at the cleavage site was further analyzed.

Detection of Gene Expression by Transcriptome Sequencing
Two days after transduction of the AdV-EGFP system and the AdV-EGFP-CRISPR system, the total RNA of CFMs and CEFs were extracted by using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), and transcriptome sequencing was performed. The concentration and quality of the extracted RNA were measured using a NanoDrop 2000 instrument (Thermo Scientific, Waltham, MA, USA) and an Agilent 2100 (Agilent, Santa Clara, CA, USA) instrument. The mRNA purification, library preparation, and sequencing steps were all performed at Shanghai Personal Biotechnology Company (Shanghai, China). A total of 12 samples (three replicates for each group) were sequenced on the Illumina Novaseq 6000 platform.
We also performed real-time PCR and chose the β-actin gene as reference housekeeping gene. ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) was used to perform the qPCR reactions in a Bio-rad CFX Connect Real-Time PCR System, with a 20 µL reaction system comprising 10 µL of ChamQ Universal SYBR qPCR Master Mix (2×), 0.4 µL of each of the forward and reverse primers (10 µM), 2 µL of cDNA and 7.2 µL of distilled water. The primers for MSTN were 1F: GCTTTTGGATGGGACTGGATT, and 1R: CAGGTGAGTGTGCGGGTATTT. The primers for β-actin were 2F: GAGAAATTGT-GCGTGACATCA, and 2R: CCTGAACCTCTCATTGCCA. The program was 95 • C for 30 s; followed by 40 cycles of 95 • C for 10 s, 60 • C for 30 s; and ended with a melting curve analysis. The 2 −∆∆Ct method was used to calculate the relative expression levels between the KO and WT groups. Statistically significant differences between groups were calculated by Student's t-test, and p < 0.05 was considered statistically significant.

Principal Component Analysis (PCA) and Differentially Expressed Gene (DEG) Analysis
The raw data were filtered using a quality control analysis, and clean reads were obtained by removing low quality reads. The clean reads were mapped to the Gallus gallus genome (GRCg6a) using HISAT2 (http://ccb.jhu.edu/software/hisat2/index.shtml accessed on 21 May 2021). We used DEGseq2 R package to perform PCA on each sample based on the expression level to identify underlying stratification across samples. Differentially expressed genes (DEGs) were analyzed using the same R package. Genes with a p-value < 0.05 and fold-change ≥2 or ≤0.5 were considered differentially expressed.

Enrichment Analysis of Differentially Expressed Genes (DEGs)
Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the DEGs were performed using the DAVID 6.8 Functional Annotation Tool (https://david.ncifcrf.gov/ accessed on 27 June 2021). Ingenuity Pathway Analysis (IPA) (Ingenuity Systems, Qiagen, California) software was employed to identify canonical pathways affected by DEGs and to determine the molecules and network interactions with MSTN by uploading DEGs and MSTN. In all analyses, the p-value < 0.05 was considered significantly different.

Generation of CRISPR/Cas9-Mediated MSTN Knockout in Chicken Cells
Immunofluorescence studies showed that the chicken fetal myoblasts (CFMs) we differentiated from primary cultures had a purity greater than 95% ( Figure 1A). We designed single guide RNA (sgRNA) sequences targeting exon 1 of MSTN, designated as sgRNA1: TGATCAGTATGATGTCCAGA, and sgRNA2: TGTGATAATCGTCTCGGTTG; these two sgRNAs are separated by about 50bp in genome. We delivered the AdV-EGFP system and the AdV-EGFP-CRISPR system into the CFMs and into chicken embryonic fibroblasts (CEFs). After 48 hours of culture, the AdV-EGGP system produces wild-type (WT) cells, while the AdV-EGFP-CRISPR system produces knockout (KO) cells. These cells were divided into four groups: CFM_WT, CFM_KO, CEF_WT, and CEF_KO ( Figure 1B-D).
Agarose gel electrophoresis of the WT cells demonstrated a 405-bp band for the MSTN PCR product, whereas the KO cells had two bands at 405 and about 351 bp ( Figure 1E), indicating effective destruction of the targeted sequence of the MSTN gene in the KO cells. Detailed analysis of the mutations produced by CRISPR/Cas9 targeting showed that 75-80% of the sequence was mutated in both types of KO cells ( Figure 1F). A strong preference (~53%) for perfect ligation of the cut sites was evident after CRISPR treatment ( Figure 1G).

Differentially Expressed Genes (DEGs) between KO and WT CFMs
Transcriptome sequencing was used to detect the mRNA expression profiles of CFM_KO and CFM_WT. Each group contained 3 samples, and 6 libraries were constructed. An average of 31 million reads was generated for all the libraries, and more than 91% of the clean reads were mapped to the galGal6 chicken genome, suggesting good sequence quality. Principal component analysis (PCA) showed clear separation between CFM_KO and CFM_WT and concordance among group replicates (Figure 2A), indicating that the similarities between the samples in each group were high, and the two groups were completely separate. Comparison of the CFM_KO and CFM_WT muscle cells revealed 296 significantly differentially expressed genes (DEGs), including 99 upregulated and 197 downregulated genes in the CFM_KO samples ( Figure 2B, Table S2). The top ten DEGs are listed in Table 1. The expression of MSTN was downregulated in the CFM KO cells compared to the WT, and we also confirmed this result with real-time PCR ( Figure 2C). Ingenuity Pathway Analysis (IPA) revealed an interaction between follistatin like 3 (FSTL3) and cyclin dependent kinase inhibitor 1A (CDKN1A) with MSTN ( Figure 2D). Furthermore, for 25 DEGs related to muscle development in CFMs, hierarchical clustering showed that those genes were related to skeletal muscle tissue development ( Figure 2E).

Enrichment Analysis of DEGs
Gene Ontology (GO) analysis performed on the MSTN KO cells to generate classification clusters identified 535 significantly enriched entries in the biological process category in the CFM_KO vs. CFM_WT comparison (Table S3). Most of the biological process-enriched items were significantly associated with cell proliferation and differentiation, muscle growth and development, signal transduction, apoptosis, the immune system, and energy metabolism ( Figure 3A). Several major GO terms related to muscle contraction (positive regulation of actin cytoskeleton reorganization, regulation of actin cytoskeleton reorganization, actin filament polymerization, actin filament organization, and actin filament-based process) were also identified.    (Table S4), including the MAPK signaling pathway (p = 0.03) and the PI3K/Akt pathway (p = 0.04). The canonical pathways were also identified by subjecting the DEGs to IPA (Table S5). Some muscle developmentrelated pathways, such as FGF signaling, calcium signaling, and the STAT3 pathway, were enriched. FGF signaling is known to promote myoblast proliferation [23], calcium signaling is related to skeletal muscle development, maintenance, and regeneration [24], and the STAT3 pathway is a critical mediator of myoblast proliferation [25]. IPA analysis identified an MSTN-related differential gene network in CFM_KO cells that was associated with connective tissue development and function; skeletal and muscular system development and function; and tissue development ( Figure 3B). As shown in Figure 3B, MSTN can act directly on Akt, or it can act on Akt through collagen type i or C/EBP family genes, such as collagen type II alpha 1 chain (COL2A1), collagen type XX alpha 1 chain (COL20A1), and collagen type IV alpha 2 chain (COL4A2). Through enrichment analysis, we found some pathways, biological processes and networks related to muscle development, which deepened our understanding of the function and molecular mechanism of MSTN in birds.

Comparison of Transcriptome Changes after MSTN KO in CFMs and CEFs
The CEF_WT and CEF_KO cells each contained 3 cell libraries. PCA revealed clusters of similar samples in this experiment (Figure 2A). The two types of cells were completely separated, and the KO group and WT group of each cell line were completely separated. The CEFs samples showed 510 DEGs, including 178 upregulated genes and 332 downregulated genes in the CEF_KO samples ( Figure 4A, Table S6). Comparison with the DEGs of the CFMs revealed 236 differential genes that were only expressed in muscle-related CFMs but not in embryonic CEFs ( Figure 4B, Table S7). Among the 25 DEGs related to muscle development in CFMs, only 3 genes were also found in CEFs, namely Potassium inwardly rectifying channel subfamily J member 5 (KCNJ5), interleukin 1, beta (IL1B), and dopamine receptor D1 (DRD1). The other 22 genes were unique to CFMs.
In total, 79 significantly enriched entries were identified in the biological process category in the CEF_KO vs. CEF_WT comparison (Table S8). The enrichment results for CEFs and CFMs indicated some common enrichment pathways ( Figure S2). Comparison of the results for CEFs and CFMs revealed an enrichment of DEGs associated with the processes of cell differentiation, cell growth, cell adhesion, cell activation, and MHC class II biosynthesis in the CFMs ( Figure 4C).
The results of KEGG pathway enrichment analysis revealed 35 pathways enriched with DEGs in the CEF_KO vs. CEF_WT comparison (Table S9). Cytokine-cytokine receptor interaction and NF-kappa B signaling pathway were present in the enrichment results of the two sets of DEGs. The MAPK signaling pathway was also enriched in CEF_KO cells (p = 0.09), but the PI3K/Akt pathway was not enriched in CEF_KO cells. Using IPA to identify canonical pathways (Table S10), we identified three typical pathways related to muscle development in CFM_KO cells, but the STAT3 pathway was not enriched in CEF_KO cells. In total, 79 significantly enriched entries were identified in the biological process category in the CEF_KO vs. CEF_WT comparison (Table S8). The enrichment results for CEFs and CFMs indicated some common enrichment pathways ( Figure S2). Comparison of the results for CEFs and CFMs revealed an enrichment of DEGs associated with the processes of cell differentiation, cell growth, cell adhesion, cell activation, and MHC class II biosynthesis in the CFMs ( Figure 4C).
The results of KEGG pathway enrichment analysis revealed 35 pathways enriched with DEGs in the CEF_KO vs. CEF_WT comparison (Table S9). Cytokine-cytokine receptor interaction and NF-kappa B signaling pathway were present in the enrichment results

Discussion
Muscle growth and development are usually regulated by specific core genes and signal transduction pathways [26,27]. Among the core genes, MSTN is the most important negative regulator [28]. The molecular mechanism of MSTN regulation in chicken muscle is not yet fully understood but understanding the potential regulatory mechanism of MSTN is very important in poultry meat production. Previous studies have used transcriptome analysis to explore the transcriptome profile of MSTN knockout (KO) muscle cells of mice and quails to understand the possible molecular mechanism, but few have explored the transcriptome profile of MSTN KO in chicken muscle cells [18,19]. Our aim in the present study was therefore to knock out MSTN in chicken muscle cells to allow a comparison of the transcriptome profiles of MSTN KO cells and WT cells. This approach systematically revealed DEGs, significantly enriched items, and KEGG pathways between the KO cells and WT cells, thereby providing reference data for the role of MSTN in muscle growth and development.
In the present research, transcriptome sequencing was used to analyze the DEGs between MSTN KO and WT CFMs. The comparison revealed TNNT1 and gap junction protein delta 2 (GJD2) as DEGs under MSTN control specifically in chicken muscle cells, in agreement with previous mammalian MSTN knockout studies in mice and pigs [29,30]. EXFABP has also been found in MSTN KO chicken leg muscle [31]. Our study also identified several previously unreported DEGs that interacted with MSTN in chicken cells. One was FSTL3, a natural inhibitor of members of the TGF-β family, which mainly binds to MSTN to inhibit its function [32]. Another was CDKN1A (also known as p21), an MSTN target gene [33].
We also performed a similar MTSN knock out in chicken embryonic fibroblasts (CEFs), a model of embryonic development, to allow us to distinguish DEGs associated specifically with muscle growth from DEGs associated with general embryonic development. We found 236 DEGs that were only expressed in CFMs but not in CEFs. Of particular relevance was the finding that 22 of the 25 DEGs related to muscle development showed differential expression in CFMs but not in CEFs. Of those 22 muscle-related DEGs, TNNT1 and EXFABP were the top DEGs in CFMs. Troponin T (TnT) is a core participant in the calcium regulation of actin filaments thin filament function, and TNNT1 encodes slow skeletal muscle TnT [34]. EXFABP is a stress protein expressed during the differentiation of myoblasts [35] and may be essential in the differentiation of multinucleated myotubes [36]. Myosin binding protein C, slow type (MYBPC1) is associated with muscle structure [37], and the ACE gene controls skeletal muscle growth [38].
In mammals, MSTN generally functions through three signaling pathways: Smad activation, MAPK activation, and inhibition of Akt signaling [13]. The MAPK signaling pathway is a major regulator of skeletal muscle development [39], whereas the PI3K/AKT signaling pathway plays a key role in promoting cell proliferation and inhibiting cell death [40]. Knockout of MSTN in the muscle-related CFMs resulted in a significant enrichment of both the MAPK signaling pathway (FGF16, CACNB4, FGF14, FGF19, CACNG5, IL1B, CACNA2D3, FGFR3, and FGF22), and the PI3K-Akt signaling pathway (CSF1R, COL2A1, SYK, COL4A2, ITGA4, HGF, IL2RG, and TLR4). By contrast, in the CEFs, knockout of MSTN resulted only in enrichment of the MAPK pathway (p = 0.09), but not the PI3-Akt pathway. This enrichment result was consistent with the results previously reported for MSTN knockout in pig and mouse cells [30,41].
Our subsequent IPA analysis revealed the network and location where MSTN participates. The network diagram was related to connective tissue development and function; skeletal and muscular system development and function; and tissue development. MSTN can act directly on the Akt in the PI3K/Akt pathway or indirectly through type i collagen or C/EBP. The IPA analysis also underscored the importance of the three main canonical pathways of DEGs in CFM_KO vs. CFM_WT involving muscle development-the FGF signaling pathway, calcium signaling, and STAT3 signaling pathway. Among these canonical pathways, the STAT3 pathway was not enriched in CEF_KO vs. CEF_WT, suggesting a specific involvement of the STAT3 pathway in muscle cell development. Interestingly, enrichment of the STAT3 signaling pathway was also reported in a MSTN knockout study in chicken leg muscles [31]. We speculate that this may be a unique regulatory pathway for MSTN in birds. These genes and pathways related to MSTN in this study provide a basis for molecular breeding to increase muscle production in commercial broilers.

Conclusions
In the present study, we knocked out the MSTN gene in CFMs and analyzed the DEGs between KO cells and WT cells to determine the role of MSTN specifically in muscle cells. We achieved about 50 bp deletions of MSTN in KO cells and clear downregulation of MSTN expression. Subsequent analysis revealed many genes that interact with or are affected by MSTN. Comparison of the transcriptome results of non-muscle CEFs with those of the muscle-related CFMs revealed the possible mechanism of action of MSTN in chicken muscle cells. Our study findings provide a valuable resource for understanding the biological functions and molecular mechanisms of MSTN for increasing muscle production in commercial meat chickens.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/genes13010058/s1, Figure S1: Schematic of the AdV-CRISPR system used in this study, Figure S2: Representative gene ontology (GO) enrichment terms of DEGs in the CFM_KO vs. CFM_WT groups and CEF_KO vs. CEF_WT, Table S1: Prime list, Table S2: The differentially expressed genes in CFM_KO vs. CFM_WT, Table S3: GO enrichment terms of DEGs in the CFM_KO vs. CFM_WT groups, Table S4: KEGG enrichment terms of DEGs in the CFM_KO vs. CFM_WT groups, Table S5: Canonical pathways about cellular growth, proliferation, and development in the CFM_KO vs. CFM_WT, Table S6: The differentially expressed genes in CEF_KO vs. CEF_WT, Table S7: The differential genes that were only expressed in CFM_KO vs. CFM_WT, Table S8: GO enrichment terms of DEGs in the CEF_KO vs. CEF_WT groups, Table S9: KEGG enrichment terms of DEGs in the CEF_KO vs. CEF_WT groups, Table S10: Canonical pathways about cellular growth, proliferation, and development in the CEF_KO vs. CEF_WT.

Data Availability Statement:
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: http://www. ncbi.nlm.nih.gov/bioproject/776535, PRJNA776535 accessed on 30 January 2022.